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We numerically analyse the rotation of a neutrally buoyant spheroid in a shear flow at small shear 
Reynolds number. Using direct numerical stability analysis of the coupled nonlinear particle-flow 
problem we compute the linear stability of the log-rolling orbit at small shear Reynolds number, 
Re a . As Re a —> 0 and as the box size of the system tends to infinity we find good agreement between 
the numerical results and earlier analytical predictions valid to linear order in Re a for the case of 
an unbounded shear. The numerical stability analysis indicates that there are substantial finite- 
size corrections to the analytical results obtained for the unbounded system. We also compare the 
analytical results to results of lattice-Boltzmann simulations to analyse the stability of the tumbling 
orbit at shear Reynolds numbers of order unity. Theory for an unbounded system at infinitesimal 
shear Reynolds number predicts a bifurcation of the tumbling orbit at aspect ratio A c ~ 0.137 below 
which tumbling is stable (as well as log rolling). The simulation results show a bifurcation line in 
the A-Re a plane that reaches A « 0.1275 at the smallest shear Reynolds number (Rc 0 = 1) at which 
we could simulate with the lattice-Boltzmann code, in qualitative agreement with the analytical 
results. 
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I. INTRODUCTION 


The angular motion of a neutrally buoyant spheroid in a simple shear has recently been studied extensively and 
in detail at moderately large shear Reynolds numbers, by numerical stability analysis and by computer simulations 
using the lattice-Boltznrann method m\- Ding and Aidun [T| analysed rotation in the flow-shear plane and found 
that a saddle-node bifurcation gives rise to steady states where the symmetry axis of the particle aligns with a certain 
direction in the flow-shear plane. The authors of Refs. [2H7j analysed this bifurcation in detail and found a large number 
of additional bifurcations at intermediate and large Reynolds numbers that give rise to intricate angular dynamics. 

At zero shear Reynolds number by contrast particle and fluid inertia are negligible, and the angular dynamics is 
determined by an infinite set of marginally stable periodic orbits, the so-called Jeffery orbits [8]. 

The effect of weak fluid and particle inertia on the angular motion of a neutrally buoyant spheroid in an unbounded 
shear was analysed recently using perturbation theory [5HTU]. In Refs. [TUI El an approximate angular equation of 
motion was derived for arbitrary aspect ratios of the spheroidal particle, valid to linear order in the shear Reynolds 
number. Linear stability analysis of the Jeffery orbits subject to infinitesimal inertial perturbations allowed to de¬ 
termine the linear stability of the log-rolling orbit (where the particle symmetry axis is aligned with vorticity), and 
of tumbling in the flow-shear plane: log rolling was found to be unstable for prolate spheroids and stable for oblate 
spheroids, in agreement with the results obtained by Subramanian and Koch fl5] in the slender-body limit. Refs. fTUlfTTl 
predicted that tumbling in the flow-shear plane is stable for prolate spheroids. For oblate spheroids tumbling was 
found to be stable for flat disks, otherwise unstable. 

This is a problem with a long history [TO]. An earlier attempt [14] to compute the stability of log rolling of nearly 
spherical particles at infinitesimal shear Reynolds number arrived at conclusions at variance with the results stated 
above, namely that log rolling is stable for nearly spherical prolate spheroids. Lattice-Boltzmann simulations of the 
problem at moderate shear Reynolds number did not find stable log rolling for prolate spheroids im but as pointed 
out in Ref. [5] the shear Reynolds number was not small enough to allow for a definite comparison with theoretical 
predictions that consider the effect of fluid inertia as an infinitesimal perturbation. 

This motivated us to analyse the stability of the log-rolling orbit numerically at small shear Reynolds number (Re a ) 
by discretising the coupled particle-flow problem directly. This method is precise enough at sufficiently small Re a to 
determine which theory is correct, and for which values of the Reynolds number it applies. We find that the theory of 
Refs. ITU} ITTl agrees excellently with the simulation results at infinitesimal Re a when the system size tends to infinity 
(the theory assumes that the shear is unbounded). Our numerical method allows to compute the effect of confinement, 
and to estimate the importance of higher-order Re a -corrections to the analytical results for the log-rolling orbit. To 
analyse the bifurcations of the tumbling orbit at small shear Reynolds numbers we use lattice-Boltzmann simulations. 
At the smallest Re a attained with the lattice-Boltzmann code (Re a = 1) the bifurcation occurs at a critical aspect ratio 
of A c ss 0.1275 in the finite system, in qualitative agreement with the analytical results obtained for an unbounded 
system. 

We briefly comment on the wider context of this paper. Recently there has been a surge of interest in describing 
the tumbling of small non-spherical particles in turbulent [T5U2T1 and complex flows [22l - f24] using Jeffery’s equation. 
Studies of the dynamics of larger non-spherical particles in turbulence Dunums] take into account particle inertia 
but neglect fluid inertia because it is difficult to solve the coupled particle-flow problem. For heavy particles this may 
be a good approximation, but the results summarised in this paper (and the results of Refs. [T» 6,IU ITUUTU]) show that 
this is approximation is likely to fail for neutrally buoyant and nearly neutrally buoyant particles. 

The remainder of this paper is organised as follows. Section |TT] describes the coupled particle-flow problem that is 
the subject of this paper. In Section [TTT1 we summarise the analytical results of Refs. ITUUTUI and find the bifurcations 
of the angular equation of motion obtained in these references. Our numerical results are described in Section |IV[ 
and compared to the analytical results. Section [V] contains our conclusions. 


II. FORMULATION OF THE PROBLEM 


The problem has the following dimensionless parameters. The shape of the spheroid is determined by the shape 
factor A defined as A = (A 2 — 1)/(A 2 + 1) where A is aspect ratio of spheroid, A = a/6 for prolate spheroids, a is 
the major semi-axis length of the particle, and b is the minor semi-axis length. For oblate spheroids the aspect ratio 
is defined as A = 6/a. The effect of fluid inertia is measured by the shear Reynolds number Re a = a 2 s/v where v 
is the kinematic viscosity of the fluid and s is the shear rate. Particle inertia is measured by the Stokes number, 
St = (p p /pt)Re a where p p and pf are particle and fluid mass densities. The numerical computations described in this 
paper are performed in a finite system of linear size L, and k = 2 a/L is a dimensionless measure of the system size, 
2 a is the length of the major axis of the particle. 
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FIG. 1. (Colour online). Schematic illustration of spheroid in a simple shear in a coordinate system that translates with 
centre of mass of the particle. Vorticity points along the negative e 3 -axis, and ei is the flow direction. The flow-shear plane is 
spanned by ei and eh- We use two different coordinate systems to express the orientation of the unit vector n aligned with the 
symmetry axis of the particle, a Spherical coordinate system used for analysing linear stability of tumbling in the flow-shear 
plane, 6 is the polar angle from the vorticity axis, and tj> is the azimuthal angle in the flow-shear plane, b Spherical coordinate 
system used for analysing linear stability of log rolling, n = [0, 0,1] corresponds to \ = ip = 0. 


We use dimensionless variables to formulate the problem. The length scale is taken to be the major semi-axis length 
a of the spheroid. The velocity scale is as, the pressure scale is /is, and force and torque scales are /isa 2 and /isa 3 
respectively, /i is the dynamic viscosity of the fluid. In dimensionless variables the angular equations of motion read 

n = wAn, St L = St (Id; + iu>) = T . (1) 

Here n is the unit vector along the particle symmetry axis, dots denote time derivatives, I = A 1 {1 — Pj_) + B 1 P_l 
is the particle-inertia matrix, Pj_ is a projector onto the plane perpendicular to n with elements Pij = Sij — nijij, 
and A 1 and B 1 are moments of inertia along and orthogonal to n. The particle angular velocity is us, and T is the 
hydrodynamic torque: 


T = 


r A tads. 


/y 


( 2 ) 


The integral is over the particle surface 5?, r is the position vector, and u is the stress tensor with elements <Jij = 
— pSij + 2Sij where p is pressure, and S tl are the elements of the strain-rate matrix S, the symmetric part of the 
matrix A of fluid-velocity gradients with elements A (? = djUi (iq are the components of the fluid velocity u). The 
anti-symmetric part of A is denoted by O with elements 0%j. To determine the torque it is necessary to solve the 
Navier-Stokes equations for the incompressible fluid: 


Re a (c^u + [u • V)tt) = — Vp + Au , V-it = 0. (3) 

For a neutrally buoyant particle Re a = St. 

It is assumed that the slip velocity vanishes on the particle surface J^, u = u) A r when The perturbation 

calculations in Refs. ITOl - fTi?! apply to a simple shear in an unbounded system, and it is assumed that the fluid velocity 
far from the particle is unaffected by its presence: u = u°° as |r| — > oo. Here u°° denotes the velocity field of the 
simple shear u°° = A°°r with = SnSj2 (see Fig. [Tj for an illustration of the geometry). The symmetric and 
antisymmetric parts of A°° are denoted by S°° and O 00 , respectively. 

The numerical computations described in this paper pertain to a finite system, a cube of linear size 2 kA 1 (in 
dimensionless variables). In the shear direction u\ = ifiT 1 at r 2 = ifiT 1 . In the flow and vorticity directions periodic 
boundary conditions are used. 


III. THEORY AT SMALL Re a . 

In Refs. HT)ljl2l an approximate angular equation of motion for a neutrally buoyant spheroid in an unbounded shear 
flow was derived, valid to linear order in Re a = St: 

n = 0 00 n+A[S 00 n-(n • S°°n)n] + /3x{n ■ §°°n)P ± S°°n (4) 

+/3 2 (n • S 00 n)0 00 n + 0 3 F ± D°° S°°n+j3 4 P ± §°° §°°n . 

The first two terms on the r.h.s. of this equation are Jeffery’s result for a neutrally buoyant spheroid in the creeping-flow 
limit. The remaining terms are corrections due to particle and fluid inertia. The four coefficients /3 a (for a = 1,..., 4) 
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TABLE I. Asymptotic behaviour of the functions b a ( A) = j3 a / Re a where f} a are the coefficients in Eq. @ for St = Re a . The 
asymptotes are found by expanding the solutions from Refs, [lp; 1111 


prolate 

oblate 

A —> 00 

A —> 0 


_7_ —197 log 2A+92 log A log 4A+106+92(Iog 2) 2 

15(2 log A —3+log 4) ' 15A 2 (2 log A—3+log 4) 2 

1 i (log A— 1+log 2)(S log 2A —7) 


5(2 log A —3+log 4) _l 5A 2 (2 log A—3+log 4) 2 


63 


4 

'5A 3 


4 

15A 2 
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30 \457T 203 ' 1357T 2 80 J 
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+ WK 


_ 1 1 ( jv_ _64_\ \ 1 / 5 _ ^1024 1 3U\ >2 

3 ^ V 20 i8ir) A '\3 135-ir 2 ' 80 ) A 


are linear in Re a and St but non-linear functions of the particle aspect ratio A: /3 a = 6 i Re “' ) (A)Re a + &i St ' ) (A)St. These 
functions were computed by Einarsson et al. pUJ, TI for general values of A, and in Ref. 1121 in the nearly-spherical 
limit. Eq. Q determines the effect of small inertial perturbations on the Jeffery orbits. It turns out that log rolling 
(n aligned with the vorticity axis) and tumbling in the flow-shear plane survive small inertial perturbations. In the 
following two Sections we discuss the linear stabilities of these two orbits, for St = Re a . We write /3 a = Re a 6 a (A). 
Table [I] gives the asymptotes of these functions for large and small values of the aspect ratio A. The asymptotes are 
obtained by expanding the results derived in Ref. ITT1 


A. Linear stability analysis of log rolling 

To analyse the stability of the log-rolling orbit we use the coordinate system shown in Fig. The angles \ an d 
ip are defined so that 


n\ = sin ip , ri 2 = cos %p sin y , 773 = cos %p cos y . 


(5) 


In these coordinates the equation of motion @ takes the form: 

=\ [4(A cos 2 ip + 1) sec ip sin y 
8 L 

+ (4/3i cos 2 ip sin 2 y + ( -2f3 2 - @3 + ^ 4 ) cos 2y + 2/3 2 + 3/? 3 + ^ 4 ) sin ip] cos ip , 

X = ^[ 2 ( A “ 1 ) tan ip + ((/3 2 - /3i) cos 2ip + /3i - fa ~ fc + Pi) siny] cosy. 

Log rolling along the vorticity direction n = [0, 0,1] corresponds to y = ip = 0, and this is a fixed point of Eq. ([6]) 
since ip = y = 0 in this direction. The stability of this fixed point is determined by the eigenvalues of the linearisation 
of Eq. ([6]) around this fixed point. To linear order in Re a the eigenvalues take the form 


( 6 a) 

( 6 b) 


Alr = ^ ± A 2 + o(Re a ). 


(7) 


The real part of this expression was derived in Refs. [10:, [IT] (see for example Fig. 3a in Ref. hod. The coefficient is 
linear in Re 0 and its sign determines the stability of the log-rolling orbit at infinitesimal Re a . The coefficient is positive 
for prolate spheroids (unstable log rolling) and negative for oblate spheroids (stable log rolling). The imaginary part 
in Eq. |7) shows that the log-rolling fixed point is a spiral at small Re a . The imaginary part has no correction to 
linear order in Re 0 . 


B. Tumbling in the flow-shear plane 


Under which circumstances is tumbling in the shear plane stable? In this Section we first summarise the results 
of analytical linear-stability calculations of Refs. IT0I - IT21 at infinitesimal Re a . Second we discuss finite but small shear 
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Reynolds numbers. To analyse tumbling in the flow-shear plane we use the coordinates employed in Refs. 11014121 
(illustrated in Fig. 


n\ = sin 9 cos (f> , n2 = sin 9 sin <j) , 713 = cos 9 . 

In these coordinates the equation of motion ([4]) takes the form 

(j) = ^ (A cos 2<f> — 1) + sin 2 9 sin Acf) — i sin 2 cf> (/^ sin 2 9 + /? 3 ) , 

2 8 4 

9 = A sin 9 cos 9 sin cf> cos ^ sin 9 cos 9 (/ 3 \ sin 2 9 sin 2 2 cf> + fj 3 cos 2 <f> + ( 3 4 ) . 


( 8 ) 

(9a) 

(9b) 


This is Ecp (42) in Ref. US 
for (j) in this plane is 


Eq. (9bI shows that 9 = 0 at 9 = 7 r/ 2 , in the flow-shear plane. The equation of motion 


4> = y (A cos 2<j> - 1) + 4/?i sin 4</> - \ (/3 2 + /3 3 ) sin 2(f>. 
2 o 4 


( 10 ) 


At infinitesimal values of Re a there is a periodic tumbling orbit in the flow-shear plane because <j> < 0. Its linear 
stability exponent yx at infinitesimal shear Reynolds numbers was calculated in Refs. US EED It was found that 
tumbling in the flow-shear plane is stable for prolate particles in this limit, and unstable for not too thin oblate 
particles. For thin platelets tumbling was found to be stable. For infinitesimal shear Reynolds numbers the bifurcation 
occurs at the critical aspect ratio [Kf; ITT] 

A c « 0.137. (11) 


This concludes our summary of the results of Refs. 1T01 [TT1 valid at infinitesimal Re a . 

As Re a increases we see that cf> > 0 in Eq. ([9J) for some value(s) of <f> . This implies the existence of fixed points in 
the flow-shear plane. This happens in Eq. © for any aspect ratio. But Eq. ([9]) is valid only to linear order in Re a . 
For this reason we only look at limiting cases where Eq. © exhibits bifurcations at small values of Re a . This occurs 
for thin rods and plates, as will be seen below. 

Consider first rods. Rods of infinite aspect ratio align with the flow direction, particles with finite aspect ratio 
tumble at infinitesimal Re a . At finite values of Re a a bifurcation may cause a rod with finite aspect ratio to align. 
To find this bifurcation point we expand cf> to second order in 1/A (Table [i]) and to second order in <f>. A double root 
of the resulting quadratic equation for (f> determines the bifurcation point: 

Re*/ 1 -* ~ ( — 3 + log 4 + 2 log A) as A —> 00 . (12) 

A 

The leading terms of this result for Re^ 1 * agree with Eq. (3.31) in Ref. |T3j (up to a factor of 87t). Subramanian and 
Koch [13] derived their result using the slender-body approximation. Note that the qualitative features of the dynamics 
in the vicinity of Re^ cl * are consistent with Eq. (12) in Ref. JT] (see also Zettner and Yoda [26]). As e = Re a — Re), cl * 
tends to zero from below the period of the tumbling orbit tends to infinity as (—e)^ 1 / 2 . Above the transition there 
are two fixed points, a saddle point and a stable node. It follows that the particle aligns at the angle 


<t> 0 


1 s/l V30 

A 15 yj A(—3 + log 4 + 2 log A) 


(13) 


for small values of e. The form of this equation is consistent with Eqs. (3.30) and (3.31) in Ref. [T5I 

Now we turn to thin disks. The symmetry vector of an infinitely thin disk aligns with the shear direction, <f> = 0 for 
< f> = 7t/2 when A = 0. For non-zero values of A the vector n tumbles in the flow-shear plane in the limit of Re a -+ 0. At 
finite (but small) values of Re a a bifurcation may cause the disk to align. To find this bifurcation point we expand 4> 
to second order in A (Table [I]) and to second order in 8<j) = (j) — 7 t/ 2 . As above a double root of the resulting quadratic 
equation for 6(f) determines the critical shear Reynolds number: 

Rei c2) ~ 15A as A -+ 0. (14) 


For Re a > Re)/ 2 * the symmetry vector n of the disk aligns in the flow-shear plane at the angle 


cf>Q = 77 + A + x“r V / 30A as A —> 0. 
2 15 


(15) 


In deriving this expression only the lowest orders in A and e were kept. 

The bifurcation lines in the A-Re a -plane given by Eqs. ( |TT] ), ( 12 ), and (14) are shown in Fig. [4] This figure also 
contains the results of our direct numerical simulations (DNS) which we discuss next. 
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FIG. 2. (Colour online), a Comparison b etween the analytical result 0 for 5R7 lr (solid red line) and numerical results from 
direct numerical stability analysis (Section IV A i. Parameters: Re a = 2.5 • 1CU 4 , k = 0.025 (circles, o), k = 0.05 (triangles, v)> 
k = 0.1 (diamonds, o), and k = 0.2 (squares, □). b Same comparison for the imaginary part Q7l R . The inset shows numerical 
results for A 7 RR for slender prolate spheroids. Shown are results for k = 0.2, Re a = 2.5 • 10~ 4 , and for different grid sizes in 
the vicinity of the particle: same resolution as in the main plot (squares, □), characteristic lengths of the finite elements close 
to the particle larger by a factor of 1.25 (<), 1.5 (A), and 2 (>). 


IV. NUMERICAL COMPUTATIONS 


We performed different types of DNS of Eqs. 0 to 0 in a finite domain with velocity boundary conditions in the 
shear direction, periodic boundary conditions in the other directions, and no-slip boundary conditions on the particle 
surface. We directly computed the linear stability of the log-rolling orbit using version 4.4 of the commercial software 
package Comsol Multiphysics™. As explained below this method could not be used to numerically determine the 
linear stability of tumbling in the flow-shear plane. Therefore we used lattice-Boltzmann simulations of the particle 
dynamics to determine the bifurcations of this orbit. To check the accuracy of the lattice-Boltzmann simulations we 
also performed steady-state DNS using version 9.06 of the commercial software package STAR-CCM+™. 


A. Direct numerical stability analysis of log rolling at finite values of Re a 

The eigenvalue solver in version 4.4 of the commercial finite-element software package Comsol Multiphysics™ m 
makes it possible to analyse the stability of the log-rolling orbit as described in this Section |28j . The symmetries of 
the problem ensure that log rolling exists not only at infinitesimal but also at finite shear Reynolds numbers. 

To determine the linear stability of this orbit it is sufficient to account for small deviations of n from the log¬ 
rolling direction n = [0,0,1], and for the fact that the particle spins around its symmetry axis. Thus we avoid 
computationally expensive re-meshing around the particle. 

The analysis proceeds in two steps. The first step is to find the steady-state solution of Eqs. ([l]) to ([3]) for a given 
value of Re a , keeping n fixed at n = [0, 0,1]. This determines the angular velocity u at which the particle spins around 
its symmetry axis. The second step is to allow for infinitesimal deviations of n and lj from this steady state. We use 
a so-called ‘arbitrary Lagrangian-Eulerian method 'M for grid refinement (deformation) close to the particle surface, 
linearise the resulting dynamics, and determine the eigenvalues of the linearised problem using the eigenvalue solver 
in Comsol, which is based on ARPACK FORTRAN routines for large eigenvalue problems Bung. The eigenvalue 
solver provides N eigenvalues 71 ,... ,7 n closest to the origin in the complex plane, ordered by ascending real parts, 
3?7i > ... > Sftq jv- 

When the shear Reynolds number is small we usually find that N— 2 eigenvalues 73 ,..., 7 jv are real (within numerical 
accuracy) with negative real parts. These are fluid modes that decay rapidly as the steady state is approached. In 
addition there is one leading pair of complex conjugate eigenvalues 71,2 with largest real part. This complex pair 
corresponds to the linear stability exponent 7 l R of the log-rolling orbit. It can have positive or negative real part, and 
the imaginary part determines the angular velocity of the particle. We must choose N large enough to ensure that 
this pair is among the N eigenvalues the solver finds. In most cases we find N = 200 to be sufficient. At larger values 
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FIG. 3. (Colour online), a Shows K 7 LR as a function of Re a for re = 0.025 and for four different values of A. The thin solid 
lines show the limiting behaviour as Re a —> 0. The thick solid lines show fits to Eq. (161. The coefficients are given in Table [TT| 
b Finite-size corrections to Kyut/Rea for Re a = 2.5 • 10~ 4 , and for the same values of A as in panel a. The thick solid lines 
show quadratic fits to the small-re behaviour, thin solid lines show the corresponding linear re-dependence for small values of re. 


of Re a it may happen that fluid modes have real parts that are larger than that of 7 RR , yet they are still real (within 
numerical accuracy). When this happens we verify that the complex pair describes the stability of the orientational 
dynamics of the particle by numerically integrating the dynamics near the steady state. 

In this way we determine 7 LR as a function of the particle aspect ratio A for different degrees re of confinement, 
and for different values of Re a . Fig. [ 2 ] shows real and imaginary parts of 7^ R as functions of the aspect ratio of the 
particle, for a small shear Reynolds number (Re a = 2.5 • 10 -4 ) and for different system sizes, re = 0.025, 0.05, 0.1, and 
0.2. Fig. [2^ compares the numerical results for the real part of 7 lr with the theory, Eq. 0. We observe excellent 
agreement for the largest system (re. = 0.025). This lends support to the analytical results of Refs. fTOMI^l and also 
to the numerical linear stability analysis. As we reduce the system size (re = 0.05, 0.1, and 0.2) we observe increasing 
deviations from the theory for the unbounded system, as expected. For re = 0.2 there are substantial finite-size 
corrections. Fig. |2ja compares numerical results for the imaginary part 37 l R with Eq. (|tJ) . Also for the imaginary 
part good agreement between the numerical results and Eq. ([7| is observed for large system sizes, at least for moderate 
aspect ratios, 10 _1 < A < 10. As for the real part there are finite-size corrections, but they are small relative to the 

0(Re°)-term in Eq. 13- 

Now consider the deviations between the numerical results and theory that can be seen in Fig. [2 Jd for more extreme 
aspect ratios. In this panel (and also in Fig. W) the size of the finite elements close to the particle surface is chosen as 
small as possible given the limited computational memory. But for very large (and also for very small) aspect ratios 
the resolution is insufficient. This can be seen in the inset of Fig. The inset shows data for for re = 0.2, and 

for different grid resolutions in the vicinity of the particle. For moderate aspect ratios the results converge quickly 
as the mesh size is reduced. But for A > 10 we do not obtain convergence, reflecting the limitations of the numerical 
approach. 

Fig. [3]a shows finite-Re a corrections to 3?7 lr for four different values of A, for the smallest value of re at which 
we could reliably compute, re = 0.025. Theory [33 El] suggests that there are Re 4/2 -corrections to Eq. ([ 7 ]) in the 
unbounded problem (re —I 0). These corrections arise as follows. The leading-order inertial perturbation of the angular 
dynamics (linear in Re a ) is obtained in terms of the solution of the lowest-order problem, Stokes problem. At finite 
but small values of Re a the Stokes solution provides an accurate description of the fluid velocity in the vicinity of 
the particle. But at larger distances from the particle (further away than the Ekman length 2a/Rey 2 ) the actual 
solutions decay more rapidly than the Stokes solution. Within the perturbative scheme used in Refs. [IOlH3j this 
gives rise to a Re^-correction. The precise form of higher-order Re a -corrections is not known. We assume that the 
next order is quadratic in Re a and compare the Re a -dependence observed in the direct numerical simulations with a 
fit of the form 


5?7lr = ai(A,re) Re a + a 2 (A, re) Re ^ 2 + a 3 (A, re) Re 2 + ... . (16) 

The values obtained for the coefficients 01 , 02 , and 03 are listed in Table |TT| The data shown in Fig. W and Table0 
are consistent with the existence of Re^ 2 -corrections when the system is large enough, re <C Re)/ 2 . 








TABLE II. Coefficients ai, 02 , and 03 from the fit of Eq. (16 1 to the data in Fig. [3^, for k = 0.025. 


A 

ffli 

a 2 

03 

1/4 

-0.0830 

0.0652 

- 0.0200 

1/2 

-0.0566 

0.0482 

-0.0183 

2 

0.0155 

-0.0125 

0.0035 

4 

0.0051 

-0.0039 

0.0008 


TABLE III. Coefficients ci, C 2 
of 64 (A)/4 to which the coefficient ci should converge as k 


and C 3 from the fit of Eq. (171 to the data in Fig. |3p. Also given are the numerical values 

0. These values are taken from Ref. ESI since the 


0 and Re a 

aspect ratios A = 1/4,1/2, 2, 4 are not small (large) enough to use the asymptotic formulae given in Table [I] 


A 

ci 64 (A)/4 

C 2 c 3 

1/4 

-0.08205 -0.0820 

0.11923 -0.04124 

1/2 

-0.05564 -0.0555 

0.09373 -0.04521 

2 

0.01526 0.0153 

-0.02784 0.01347 

4 

0.00510 0.0051 

-0.00870 0.00367 


Fig. [ 3)3 shows finite-size corrections to 3?7 lr for Re a = 2.5 • 10 4 
fits of the form 


and for four different values of A. Also shown are 


^7LR/Re a = ci(A)+ c 2 (A)k + c 3 (A)k 2 


(17) 


The resulting coefficients are given in Table a We see from Fig. that the fits describe the numerically observed 
finite-size dependence accurately, but Eq. (17) is just an ansatz. Also shown are linear approximations valid at small 


k. We see that the finite-size effects are to a good approximation linear in k for the data shown for n < 0.1. Table III 
shows that the limiting values, Ci, obtained as k —> 0 are in excellent agreement with the theoretical results for the 
unbounded system. 


B. Time-resolved lattice-Boltzmann simulations 


To analyse the bifurcations of the tumbling in the flow-shear plane we use the lattice-Boltzmann method with 
external boundary force [32] . To restrict the computational time, the domain size is set to a maximum of 240 lattice 
units. This allows us to resolve the particle with at least six fluid grid-nodes along its smallest dimension, with system 
size k = 0.2. These choices limit the range of aspect ratios that can be simulated to A € [1/8, 8]. We take Re a larger 
than or equal to unity in our simulations. This is because it is computationally very expensive to reach small values 
of the shear Reynolds number (as discussed by Rosen et al. Ena). 

To estimate the critical aspect ratio A c where tumbling changes stability for oblate particles we proceed as follows. 
We initialize the particle at rest, close to the tumbling orbit at <j> = 7t/ 2 and 0 = 7r/2 —69 with SO = 0.017. We integrate 
the dynamics for aspect ratios A = 1/8,1/7,1/6,1/5,1/4, and for Re a between 1 and 10 with unit increments. We 
determine whether the trajectory tends to tumbling in the flow-shear plane or to the log-rolling orbit, and determine 
the location of the bifurcation by interpolation. At Re a = 1 we run simulations for A ranging between 0.125 and 
0.160 with increments of 0.05 and determine the bifurcation point by linear interpolation. The results are illustrated 
in Fig. [I] We see that the results agree fairly well with Eq. (11). At the smallest value of Re a simulated with the 
lattice-Boltzmann code, the transition occurs at A c « 0.1275, not too far from the analytical result (11) at infinitesimal 
Re a for the unbounded system. 


Using lattice-Boltzmann simulations we also obtain estimates for Re)/ 1 ' 1 and Re)/ 2 -* (Section III). This is done by 


initialising the particle at rest at < j > = 7t/4 and 0 = 7r/2forA>l and at < f > = 37t/4 and 0 = 7t/2 for A < 1. We then 
determine whether the particle tends to a steady state or continued to tumble, and determine the critical Reynolds 
number by linear interpolation. The results of these simulations are also shown in Fig. [4] and are compared with the 
analytical results for thin disks and rods given by Eqs. (12) and (fl4). We find that the agreement is only qualitative. 
This is not surprising since Eqs. (12) and (141 are based on Eq. (4} that is valid only to linear order in Re a and 


cannot be expected to describe the dynamics at Reynolds numbers of order unity or larger. We also note that the 
lattice-Boltzmann simulations were performed for a rather small system, while the analytical results pertain to an 
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FIG. 4. (Colour online). Bifurcations of the tumbling orbit in the flow-shear plane. Bifurcation lines derived in Section 
III for the unbounded system - Eqs. (121, (141, and (111 - are shown as solid lines. The label TS indicates that tumbling is 


stable, TU that it is unstable, and FP that the tumbling orbit has bifurcated giving rise to a fixed orientation in the flow-shear 
plane. The dashed line denotes the symmetry line at A = 0 where the tumbling orbit changes stability. Numerical results 
for the finite system (k = 0.2) are shown as symbols: circles (o) denote results from the time-resolved lattice-Boltzmann 
simulations described in Section [IVB| crosses (x) represent results from the steady-state simulations described in Section [lV C| 
The bifurcations where tumbling in the flow-shear plane changes stability are shown in red, the bifurcations where stable 
tumbling in the flow-shear plane changes to a stable fixed point are shown in blue. 


unbounded system. Figs. [2^ and [Sfo show that there are substantial finite-size corrections to the stability exponent 
of the log-rolling orbit in the finite system, for n = 0.2. We therefore expect that there are equally important finite- 
size corrections to the locations of the bifurcations in Fig. [4] But at present we cannot perform lattice-Boltzmann 
simulations for larger systems with sufficient resolution to quantify this statement. 

In order to check the accuracy of the lattice-Boltzmann simulations at k = 0.2 we determined the critical Reynolds 
numbers Ilej !)' 11 and Re^ c2 ^ using an alternative approach, described in the next Section. 


C. Steady-state simulations using STAR-CCM+™ 

We compute the critical Reynolds numbers Re^ c1 ^ and Re^ c2 -* using version 9.06 of the commercial finite-volume 
software package STAR-CCM+™ [53], We choose the same system size as in the lattice-Boltzmann simulations, 
k = 0.2. The particle orientation is fixed at 9 = tt/2, (f> £ [0, 7t/ 2] for prolate particles, and (j> £ [7t/2,7t] for oblate 
particles. For a given particle aspect ratio A and value of Re 0 we compute the steady-state torque on the particle. 
If the torque vanishes, the chosen particle orientation is a fixed point for the given parameters. A fixed particle 
orientation makes it possible to use a very fine local grid around the particle. For different choices of <j> we find critical 
Reynolds numbers where the steady-state torque vanishes. The minimum of this critical Reynolds as a function of (j> 
gives Re^ or Re„ c2 \ for prolate and oblate particles respectively. The corresponding results for Re^ and Re^ 
are also shown in Fig. [4] We conclude that the lattice-Boltzmann simulations slightly underestimate the critical value 
Re„ cl \ while they slightly overestimate Re^ c2 \ 


V. CONCLUSIONS 

Using numerical linear stability analysis we computed the stability of the log-rolling orbit of a neutrally buoyant 
spheroid in a simple shear at small Re 0 . For infinitesimally small Re a in the unbounded system this problem was 
recently solved for arbitrary aspect ratios using perturbation theory in the shear Reynolds number. The fact that both 
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calculations agree in the limits Re a —> 0 and k —> 0 (unbounded system) lends support to the analytical calculations 
[TOHE], but also to the numerical linear stability analysis described in the present article. In the limit of large system 
size (k 0) we found that there are corrections to the analytical result for the exponent 3?7 lr that are consistent 
with terms of order Re 3 / 2 . We also investigated finite-size corrections to 3?7 lr at small Re a , and found that they are 
substantial. It would be of interest to calculate both finite-Re a and finite-size corrections to 5ft 7 lr by extending the 
method used in Refs. HMH 

We did not investigate the stability of the tumbling orbit with numerical linear stability analysis because the 
required re-meshing is computationally very expensive. Instead we studied the stability of tumbling in the flow-shear 
plane using lattice-Boltzmann simulations. We tracked the bifurcation line between stable/unstable tumbling for thin 
oblate spheroids (solid red line in Fig. [4]) down to as small values of Re a as we could reliably achieve and found that 
the transition occurs at A c ~ 0.1275 at Re a = 1, in fair agreement with the theoretical prediction 0.137. 

Finally we determined for which values of A and Re a tumbling in the flow-shear plane bifurcates to a fixed point, 
using lattice-Boltzmann simulations, and also by numerically computing steady-state torques using STAR-CCM+™. 
The two numerical procedures give results that are in fairly good agreement with each other, yet the agreement with 
the analytical results (121 and (14) is only qualitative. 

Detailed analysis of the lattice-Boltzmann dynamics near the bifurcation at Re^ 2 -* reveals the phase-space topology 
near the bifurcation at moderate Reynolds numbers (Re^ c2) ~ 7.8 at A = 1/4), see Fig. 3(d),(e) in Ref.0 For A = 1/4 


a second transition occurs at Re*/ 3 -* ss 5 where the log-rolling orbit changes from stable spiral to stable node. Eq. 
also exhibits this transition. But since Eq. Q is valid to linear order in Re a the bifurcation can only be analysed in 
the limit A —> 0. We find that the two transitions occur in reverse order: the tumbling—^fixed point bifurcation occurs 
before the spiral—mode transition as the shear Reynolds number is increased. There are several possible explanations 

41) such as the Re 3 / 2 -corrections 


for these subtle differences. They could be due to higher-order Re a -corrections to Eq. 
alluded to above. But we have also observed (not shown) in the numerical simulations of the bounded system that 
Rei c3) increases as k becomes smaller. In the limit of k —> 0 we except that the order of the transitions agrees with the 


prediction for the unbounded system. In summary we can conclude that the results of our numerical computations 
agree well with the theoretical predictions at infinitesimal Reynolds numbers: we find excellent agreement for the 
stability exponent of the log-rolling orbit, and the bifurcation of the tumbling orbit for thin oblate particles occurs 
in both theory and simulations, at similar values of A c . But there are a number of subtle differences between theory 
and simulations at larger Reynolds numbers. At present we cannot reliably perform lattice-Boltzmann simulations 
at much smaller Reynolds numbers than those shown in Fig. [4j and it is very difficult to perform such simulations 
at still smaller values of n. Therefore it would be of great interest to extend the analytical calculations to include 
Re 3/,2 -corrections, and to account for finite-size effects. 
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